function y = GAM (x, p, p0)

   ## usage:  y = GAM (x, p, p0)
   ##
   ## probit transform of gamma function

   if (any(p<=0))
      y = inf(size(x)) ;
      return ;
   endif
   y = gamcdf(x, p(1), p(2)) ;
   y = (1-p0) * y + p0 ;
   y = norminv(y) ;
   y(!isfinite(y)) = inf ;
#   y(isinf(y)) = nan ;

endfunction
